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1. Introduction. A need for a fast and accurate steady-state compressible flow 
solver exists in many areas of science and engineering. Perhaps, such a need is par- 
ticularly acute for the problem of aerodynamic design. In this case the steady-state 
solutions of the compressible flow past bodies have to be computed repeatedly, each 
time with variations in the body’s geometry. This allows to find the shape with op- 
timal in a certain sense aerodynamic parameters relying on computations only. Thus 
the necessity of the costly wind-tunnel experiments can be largely reduced. 

The search for the genuinely multidimensional schemes for the compressible Euler 
equations was motivated by the expectation, that they will provide new possibilities 
(comparing to the dimensionally split schemes, which are of the common use now), 
like: 

• capturing physics of the fluid flow more accurately; 

• constructing a more efficient steady-state (multigrid) solver. 

Considerations regarding the first point can be found in [14], [15]. The main motivation 
for constructing the truly multidimensional scheme in this work is the improvement 
of multigrid efficiency. It was observed in [20] that pointwise Gauss- Seidel relaxation 
is unstable when applied directly to the high-resolution dimensionally split scheme 
even in the simple case of two-dimensional scalar advection equation. Therefore, the 
steady-state Euler solver constructed in [21] relies on the defect correction technique, 
which is not a fully efficient way to utilize multigrid methods. Another possibility 
to avoid this difficulty is to use a multigrid algorithm that employs a well-known 
multi-stage Runge-Kutta relaxation technique which was developed in [9], [7], [8]. 

However, any further improvement of the multigrid efficiency requires to address 
this problem directly: it is necessary to develop a new high-resolution (at least at 
the steady-state) discrete scheme, such that Gauss-Seidel relaxation is stable when 
applied directly to the resulting discrete equations. This was the main motivation for 
constructing the genuinely two-dimensional advection scheme in the control volume 
context in [17], [18]. The novel feature of this scheme was the two-dimensional limiter, 
i.e. the limiter function that relies on the ratio of two finite differences in different 
directions. 

There has also been developed another class of genuinely two-dimensional ad- 
vection schemes - the so-called “fluctuation-splitting” schemes, (see [6], [22]). These 
schemes are also equipped with a nonlinear mechanism that allows to combine high- 
resolution and positivity property. Several variants of such a nonlinear mechanism 
were devised using some geometric considerations. The remarkable feature of this 
approach is the simplicity of the schemes formulated on unstructured triangular grids. 

The strong relationship between the fluctuation-splitting and control volume ad- 
vection schemes was established in [19]. As a result, the fluctuation-splitting scheme 
that utilizes the two-dimensional limiters was constructed. The action of the limiter 
function can be given a purely algebraic interpretation. 

Even though some genuinely two-dimensional advection schemes where available 
already several years ago, the task of extending these ideas to the systems of equations 
appeared to be rather difficult. 

The “algebraization” of the advection scheme formulation appeared to be crucial 
for the purpose of this work - constructing a truly multidimensional scheme for the 
Euler system. The resulting scheme is capable of producing high-resolution steady- 
state solutions. Its unique feature is that the Gauss-Seidel relaxation is stable when 
applied directly to the high resolution discrete equations. 
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The paper is organized as follows: in §2 we give a brief introduction into the 
fluctuation-splitting approach for the scalar advection equations and present consider- 
ations for extending of this approach to the systems of equations in two dimensions. In 
§3 we take a closer look on the Roe scheme for the Euler equations in one dimension. 
Then we construct a truly two-dimensional scheme for the Euler systems on structured 
(in §4) and unstructured (§5) triangular grids. In §6 we outline first as a preliminary 
the construction of the three-dimensional advection scheme. Then we present a truly 
three-dimensional scheme for the Euler system. §7 describes the multigrid algorithm 
that employs the constructed truly two-dimensional scheme. Numerical experiments 
are presented in §8, followed by discussion and conclusions in §9. 

2. Preliminaries and motivation. In the first part of this section we illustrate 
the fluctuation-splitting approach, first on the example of the scalar advection in one 
dimension. Then we present the construction of a simple truly two-dimensional advec- 
tion scheme that relies on two-dimensional limiters. In the second part we discuss the 
difficulties that arise when trying to apply the scalar advection schemes to discretize 
the systems of equations in multidimensions. 

2.1. Fluctuation-splitting approach. We shall give here a brief description 
of the fluctuation-splitting advection scheme in one and two dimensions (for details 
see [6], [19]). These schemes are required to have properties of positivity and linearity 
preserving. We shall define here the positivity property. 

Definition 2.1. A scheme is said to be of the positive type if any solution value 
on the new time level obtained by this scheme can be written as a positive combination 
of the values from the previous time level. 

Solutions obtained by the means of positive schemes will satisfy a certain maxi- 
mum principle and, therefore, do not exhibit oscillatory behavior in presence of dis- 
continuities. 

The notion of linearity preserving is used to characterize high-resolution at the 
steady-state schemes. It is trivial in one dimension. Therefore, we shall give its precise 
definition in §2.1.2. 

2.1.1. Advection in one dimension. Consider a linear one- dimensional ad- 
vection equation 

(1) u t + au x = 0 

We are interested in solving it numerically on a grid with meshsize h (see Fig.l). 
Fluctuation is defined as the residual of the equation on the linear element multiplied 
by the volume of the element. In our case the fluctuation on the segment \i — 1 ,i\ is 
defined as the residual of (1) on this segment multiplied by its length h 

R = —a(ui — Ui_i) 

The numerical solution can be interpreted as a two-stage process 

• compute the fluctuation on each element and distribute (split) it among the 
vertices of this element; 

• update the solution at each vertex due to the accumulated portions of fluctu- 
ations. 
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In our case if the fluctuation on the segment \i — 1 ,i\ is distributed equally between 
the adjacent gridpoints 

hu^l =hu n l _ x +^R + C'.O.E. 
hu? +1 =hu f +-R + C.O.E. 

we obtain the central scheme, which is neither positive nor stable, though second order 
accurate in space. Here LL C.O.E.” stands for “contributions from other elements” and 
will be omitted in the remainder of this paper. 

The upwind scheme can be obtained by adding a proper amount of the artificial 
viscosity to the central scheme. The fluctuation distribution formulae thus become 

hu£t = huU + T 1 (R-R) 
hu^ +1 = hu^ +1(R + R), 

where 

R = — \a\(ui — Ui_ i) = sign (a)R 


is the artificial viscosity. 

It is easy to see that in this case the entire fluctuation on the segment [i — 1 ,i\ 
contributes either to the solution update either at the left or at the right nodes of the 
segment depending on the advection speed direction (sign). 

2.1.2. Two-dimensional advection. Consider a linear two-dimensional advec- 
tion equation 


u t + au x + bu y = 0 

The fluctuation on the triangle T (see Fig. 2) is given by 

(3) R = R X + R\ 

where 

R x = — | [a(u 0 - u 3 )] 

Rv = [b(u 3 - u 4 )\ 

Assuming that the fluctuation distribution formulae in this case are 

h 2 u% +1 = h 2 uS + i(R x + R x ) 

(4) h 2 u™ +1 = h 2 u!i + U(R X - R x ) + (Ry + Ry ) ] 
h 2 ui +1 = h 2 ui + ^{Ry - Ry) 

where R x , R y are the artificial viscosity terms defined by 

, s R x = sign^)!?^ 

R y = sign (b)R y , 

we obtain the dimensional upwind scheme which is positive , but only first order accu- 
rate. 
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Definition 2.2. The fluctuation-splitting scheme is called linearity preserving if 
whenever the fluctuation on the triangle T vanishes then the scheme leads to a zero 
update in each of the three vertices of the triangle. 

It was observed in [6] that the linearity-preserving schemes on structured grids 
produce numerical solutions that are second-order accurate in the steady state. 
Introduce the following quantities 

R x * = R x + Ry^{Q) 

( 6 ) R y* = R y +R xmi 


where 


(7) 


Q 


R x 

Ry 


and if is a Lipschitz continuous limiter function such that 


(8) 


0 < ®(Q) <1, 0 < 


m) 

Q 


< 1 


and 


(9) *(1) = 1 

Substituting R x *,R y * for R x ,R y into (4) and (5) we obtain a linearity preserving 
(second order accurate at the steady-state in the case of structured grid) scheme. 
Using the following identity 

( 10 ) R y V(Q) = -R x ^p- 

it is easy to see that the scheme defined by (4), (5) and (6) is of the positive type. 

It is also obvious from (10) that such scheme is conservative because 

R x * +R y * =R x + R y = R 


(for more details see [19]). 

2.2. Hyperbolic systems of equations. We shall explain here what is the 
main obstacle encountered when constructing a truly two-dimensional numerical scheme 
for a hyperbolic system of equations. We shall also describe what is our approach 
to overcome it. The understanding of the basic differences between one- and two- 
dimensional cases is required for this purpose. These differences will be illustrated on 
the linear systems of equations. 

2.2.1. O ne dimension. Consider a system of conservation laws in one dimen- 
sion 

(11) u t + f(u) x = 0 

For the purpose of the discussion here it is sufficient to look at the linear constant 
coefficient case of (11). Consider a the following system of partial differential equations 


u t + Au x = 0 


( 12 ) 
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where A is N x N matrix. The system (12) is said to be hyperbolic if matrix A has a 
complete set of real eigenvalues. 

Denote T - the matrix of right eigenvectors of A. Then 

A = T _1 AT 

is a diagonal matrix. Introducing characteristic variables 

w = T -1 u 

(12) can be rewritten as a set of N decoupled advection equations: 

(13) w t + Aw x = 0 

- a set of N decoupled advection equations. It is clear from (13) that a one- dimensional 
advection scheme can be applied in a straightforward way to solve a (linear) hyperbolic 
system of equations in one dimension. A discussion concerning the nonlinear Euler 
system in one dimension will be presented in the §3. 

2.2.2. Two dimensions. A linear system of partial differential equations in two 
dimensions of the following form 

(14) Uf + Au x + Bu y = 0 
is said to be hyperbolic if the matrix 

A = cos cf)A + sin cj)B 

has a complete set of real eigenvalues for \/(f> : 0 < cj) < 27T. In this case there exist 
matrices T 4 and Tg such that Tj " 1 AT a and T^BTg are diagonal. This is usually 
utilized to construct the so-called dimensionally split schemes. 

However, matrices A and B in general do not commute. Therefore, T 4 7 ^ Tg, i.e. 
a hyperbolic system in two dimensions cannot be represented as a set of decoupled 
advection equations (unlike in one- dimensional case). This means that a truly two- 
dimensional advection scheme cannot be applied in a straightforward way to discretize 
a hyperbolic systems of equations. 

Much of the research effort in the last several years concentrated on Ending a 
way to apply the multidimensional advection schemes to discretize the hyperbolic 
systems of equations, in particular to the Euler equations of gas dynamics. One of 
the major directions was the so-called wave modeling (see [14], [5]). This approach 
concentrated on finding a way to represent (locally) the physics of two-dimensional 
flow of a compressible fluid by a finite number of simple waves, each one having an 
associated advection equation. 

The approach we present here is concerned not with applying an advection scheme 
to discretize a system of equations in two dimensions, but rather with applying to the 
systems of equations the same strategy that was used when constructing a scalar 
advection scheme. The construction of the genuinely two-dimensional scheme for 
the Euler will be presented in §4 for the structured grids and in §5 for the general 
unstructured grids case. 
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3. Euler equations in one dimension and Roe scheme. Euler equations of 
gas dynamics in one dimension can be written as follows 


u t + F(u) x = o, 


where 


u = pu 


the enthalpy H is defined by 


the speed of sound 


pu \ 

F(u) = pu 2 + p , 

V puH I 


e + p c 


b — it , 

7-12 ’ 


and the pressure 

(19) p = (7 - l)(e - ^pu 2 ) 

The Jacobian matrix 

, , A dF 

< 201 A = 

has real eigenvalues only. 

In order to apply the upwind differencing approach to the nonlinear system (15) 
a particular conservative linearization procedure was introduced by Roe [13]. For two 
states U[ and u r a matrix A(u[,u r ) was constructed having certain properties called 
together Property U. We address the reader to [13] for the details. Let us only mention 
here one of these properties: the following identity holds for any U(,u r : 

(21) A{u h u r ) ■ (u r - ui) = F(u r ) - F(ui). 

3.1. Fluctuation-splitting formulation. Define the fluctuation on the consid- 
ered element (segment \i — 1 ,i], see Fig.l) 

(22) R = ~[F(ui) - F{u t _ 1 )]. 

Following the Roe-linearization procedure we assume that the quantities which vary 
linearly on the segment [i — 1 ,i\ are the elements of the so-called parameter vector 
(see [13]) 

(23) m = (mi, ra 2 , m 3 , m 4 ) T = yfp{l, u, v, H) T . 

Therefore, its average value on \i — 1, i\ is 


m 8 _i + m, 
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Thus the following Roe-averaged quantities can be defined 

,2 

(25) 


p = TOJ 

u = m,2/mi 
H = /rhi 


and 

(26) 


c' 2 = ( 7 - 1 )(H - -) 


Introducing the Roe-averaged Jacobian A = A(v,i-i,Ui) and recalling (21) we can 
rewrite the fluctuation (22) in the following form 


(27) 


R = —A ■ (ui — Ui-i 


Having in mind the construction of an upwind scheme, it is convenient to use the 
following representation of the matrix A (see [13]) 


(28) 

where A is a diagonal matrix 

(29) 

where 


A = TAT” 1 , 


A = 


/ Ai 0 0 

0 A 2 0 

y 0 0 As 


Ai = u, A 2 ,3 = u±c 

~ \ ~ 2 ~ 3 

are the eigenvalues and T = (E , E , E ) is the matrix of right eigenvectors of A 


(30) 


E 1 = 


( 1 


E 2 = 


V “72 


/ 1 

V 


u — c 
H — uc 


e 3 = 


( 1 

u + c 
y H + uc 


The fluctuation-splitting formulation of the Roe scheme can be written as follows 

/ op hu^ = huU +UR+R) 

[ ] hu/ +1 = hu/ +1(R-R), 

where 

(32) R = -f\A\(w t - m 8 _i) 
is the artificial viscosity, and 

(33) w = T~ x u 

are the so-called characteristic variables. 

Alternatively, the artificial viscosity can be also written as 


(34) 

or 


R = -\A\{u t - Ui_i) 


(35) 


R = -sign(7)i? 
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3.2. Roe scheme revisited. It is convenient for the purpose of this work to 
introduce the following auxiliary variables (s,n,p) T , where 

(36) ds = dp ^ 

c 2 

is the entropy. The non- conservative formulation of the Euler equations in these 
variables takes the following form 

St + US X = 0 

(37) pu t + puu x + p x = 0 

Pt + Wx + pc 2 u x = 0 

Introduce the fluctuations of the Euler equations in the auxiliary variables 



( 45 ) 



or 


(46) r = sign(A)r 

It is evident from (43) and (45) that the matrix of right eigenvectors of A is 

f = C~ x f 


or T = (ei,e 2 ,e 3 ) where 


(47) 


e i = 


( 1 

0 

V° 


( 0 


e 2 = 


( 0 


— c 
~2 


e 3 = 


C 

~2 


Remark 3.1. It is obvious from the structure of the matrix T that the auxiliary 
variable s coincides with the first element of the characteristic variable vector (33) 


(48) 

Denote 


Wi = s 


M = sign(A) 


In order to make a computer program more efficient it is useful to write an explicit 
expression for the matrix M . It is easy to see that in the supersonic case 

(49) M sup = sign(h) • / 

where / is a 3 x 3 unity matrix. 

Some algebra reveals that in the subsonic case matrix sign(A) has a very simple 
structure as well 


(50) 


M sub 


t sign(h) 0 0 \ 

0 0 1/7 

v 0 7 0/ 


The second and third rows mean that the artificial diffusion added to the momentum 
equations is proportional to the fluctuation of the pressure equation and vice versa. 

Remark 3.2. The auxiliary variables formulation (37) of the Euler equations 
is, perhaps, the simplest way to write the Euler equations. The usefulness of this 
formulation will become even more evident in §§^,h. 

4. Construction of the two-dimensional Euler scheme. Euler equations of 
gas dynamics in two dimension can be written as follows 


(51) 

where 


( 52 ) 


( P \ 

pu 

pv 

V e / 


u t + F(u) x + G(u) y = o, 



( pu \ 


f pv ^ 

F(u ) = 

pu 2 + p 

■ G(u) = 

puv 


puv 


pv + p 


\ puH ) 

V pvH ) 
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where the enthalpy H is defined by 


(53) 

the speed of sound 

(54) 

and the pressure 

(55) 


e + p 
P 


7 - 1 


+ 


u 2 + v 2 



p = (n - !)( e - p 


2 , 2 
W + V s 


The quasilinear non- conservative formulation of the Euler system in auxiliary variables 
( s,u,v,p ) can be introduced in two dimensions as well 


(56) 


St + US X + vs y = 0 

put + puu x + pvUy + p x = 0 

pv t + PUV X + pVVy + Py = 0 

Pt + Wx + VPy + pc 2 (u x + Vy) = 0 


where ds = dp — . 

Remark 4.1. Note that the entropy (s ) evolution is subject to the two-dimensional 
advection equation, which is locally decoupled from the rest of the system. 

The fluctuation of the system (51) defined over the triangle T 


(57) 


R = 


u t = - 


F x + Gy) dx dy = -S T F x + G y 


where F x ,G y are some averaged values of the flux derivatives over the triangle T. 

Our construction of the truly two-dimensional Euler scheme utilizes the two di- 
mensional extension (Roe, Struijs and Deconinck [16]) of the Roe conservative lin- 
earization for one- dimensional case. Therefore, following the procedure developed in 
[16], we assume that the quantity which varies linearly over an element is the “param- 
eter vector” 


(58) m = yPp(l,u,v,H) T 

and its averaged value on the triangle T is given by the following 


(59) 


m = 


mi + m> + m 3 

3 


Roe-averaged quantities can be introduced 



u = 

(60) 

v = m 3 /rhi 
H = m^/rhi 

and 

(61) 

c 2 = (7 - 1 )[H ~ i( u- 
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Fluctuations of the Euler system in the auxiliary variables can be presented as 


(62) 


r = r x + r y , 


where 

r x = -S T A ■ {^,]m x ,pF x ,pl) T 

r y = -S t B ■ (^^ y ,pVy,p^) T 

with 


{ 

u 

0 

0 

° ^ 


/ 

V 

0 

0 

0 \ 


0 

u 

0 

1 

, B = 


0 

V 

0 

0 


0 

0 

u 

0 


0 

0 

V 

1 

\ 

0 

7 2 

0 

u ) 


V 

0 

0 

c 2 

7 j 


and St = h 2 / 2 is the area of the triangle T, and 


(63) % 

(64) pvT x 

(65) pvi 

(66) 


2m 1 (m 1 ) x 

mi(m 2 ) x - m 2 {mi) x 
mi(m 3 ) x - m 3 {mi) x 

7 — 1 

[(m 4 (mi)i + mi(m 4 )i) + (m 2 (m 2 ) x + fn 3 (m 3 ) x )} 


The corresponding terms involving derivatives in y direction can be written in the 
analogous manner. 

Introducing the matrix 


(67) 


C a 


/ 1 0 0 
u 1 0 

7 0 1 

y ( u 2 + 7 2 ) /2 U V 


1 /c 2 \ 

u/c 2 
7/c 2 

l/( 7 -l) + (7 2 + 7 2 )/(27 2 ) / 


we can define 


(68) 


iZ* = C a r x 
R y = 0" 


It can be easily verified that 
(69) 


R x = -S T JAc 
Ry = - s T Gy , 


where F x ,G y are the same averaged flux derivative values as defined in [16]. It is also 
obvious that the entire fluctuation 


(70) R = R x + R y = C a (r x + r y ) = C a r 

Consider triangle T as illustrated on Fig. 2. The fluctuation is distributed accord- 
ing to the following formulae: 



Su^ +1 

= Sul 

+ ^C a (r x -r x ) 

(71) 

Su^ +1 

= Sul 

+ T 2 C a [(r x + f x ) + (r 


Su^ +1 

= Sul 

+ ^C a (r y + f y ) 
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Defining 


r x = sign(A)r x 

[ ) f y = sign(5)r*« 

we obtain the scheme that is similar to the standard Roe dimensionally split scheme. 
The only difference is in the linearization procedure. 

We can construct now a (linearity preserving) second order accurate scheme. First, 
we shall introduce vectors r x * ,r y * with their elements defined by 


(73) 

for i = 1,2, 3, 4, where 

(74) 


r f = r f + 

r y* _ r y , 

i i ' n: i 


ft = 


and T is a (non-compressive) limiter. 

Substituting r x ,r y for r x ,r y in (71) and (72) we obtain a linearity preserving 
(second order accurate in this case) scheme. This is because if 


r = 0 


then 


r x * =r y * = 0 

as well. Therefore, the update of any variable in all the three vertices of the triangle 
T due to the fluctuation on this triangle is zero. 

The resulting scheme is also conservative since 

r = r x + r y = r x * + r y * 


Denote 


It is interesting to note that 


M x = signM) 
M y = sign (B) 


(75) 


M s x u \ if \u\ < c 

M x up , if \u\ > c, 


and 

(76) 


My 


My Ub , if |7| < c 

My Up , if |7| > c, 


where 


(77) 


t sign(h) 0 0 0 ^ 

0001/7 
0 0 sign(h) 0 

V 0 7 0 0 / 
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/ sign(7) 0 0 0 \ 


(78) 

My Ub = 

0 

0 

sign(7) 0 

0 0 

0 

1/7 ’ 


V 

0 

0 7 

o 

and 





(79) 


M s x uv 

= sign(h)/, 


(80) 


M s y uv 

= sign (6)/. 


5. 

General unstructured 

grid formulation. 

Consider an unstructured tri 


angular grid covering out domain and triangle T belonging to this grid (see Fig. 4). 
Denote hi length of the ith face of this triangle, 7) a unit vector along the ith face 
in clockwise direction. Denote also n 4 - a unit inward normal to the ith face and St - 
area of the triangle T. In this section we shall describe first the general (unstructured 
grid) version of fluctuation-splitting advection scheme. Then we shall construct the 
genuinely two-dimensional numerical scheme for the compressible Euler system, also 
formulated for the unstructured meshes. 

5.1. Advection scheme. Consider a linear constant coefficient advection equa- 
tion 


(81) u t + A • Vti = 0 

Consider a non-orthogonal coordinate system £, ?/ aligned with the vectors ei,e 2 (see 
Fig. 4) and introduce the following quantities 


(82) 

where 

(83) 


a = 


A • n 2 


P 


A • h\ 
9 


9 = e*i • n 2 = e 2 • rai . 


Then Eq.(81) can be written in the coordinate system £,7/ 


(84) ut + an?: + [3u, q = 0 
The fluctuation of the Eq.(81) 

(85) R = R ( + R\ 


where 


/o fi \ Rt = -S T [a(u 3 - ni)/hi] 

1 J BP = -St[{3(u 2 - u 3 )/h 2 ] 

The fluctuation distribution formulae 

s lU ^ +1 = + !(£« - M) 

(87) S 2 u% +1 = S 2 u% + %(BP + BP) 

S 3 ul +1 = S 3 u^ + f [(# + Rt) + (RV- BP)] 
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with the artificial viscosity terms define as follows 


(88) 

= sign (a)R^ 
R' n = sign (fi)R v 

result in a positive type upwind scheme. Introduce the following quantity 

(89) 


Define 


(90) 

B? = Rt + K n V(Q) 

RV* = RV + Ri^W 

and substituting R^* ,R V * instead of R^,R V in (87), (88). Using the identity 

(91) 

M(Q) = 

we can rewrite R^* ,R V * in 

the following form 

(92) 

R c = R c (1 IM) 
RP* = RP{ 1 - T(Q)) 

it is easy to verify that the 

resulting scheme is positive if 

(93) 

0 < ®(Q) <1, o < -5-^ < 1 

and linearity preserving if 
(94) 

Note also that 

T(l) = l 

(95) 

R = R e +R ' n * = R* + R' n . 


Therefore, the constructed scheme is conservative. 

There are two differences between the formulated scheme and the typical fluctu- 
ation splitting advection schemes (see [6], [22], [19]) 

1. no distinction is made between triangles with two outflow faces (Type I) and 
one outflow face (Type II) 

For the advection problem this implies more computational work since limiter 
function has to be evaluated on all the triangles (not only on those of Type 
II). However, for the Euler system no significant savings can be made by this 
distinction, since limiters have to be used in on all the triangles anyway (see 
§5.2). 

2. the derivatives are approximated by finite differences along the faces which 
were chosen arbitrarily and are not necessarily both inflow or outflow (a and 
(3 are not necessarily of the same sign). 

In case of the advection equation choice of the faces which are both either 
inflow or outflow for derivative approximation provides better resolution of 
discontinuities. For the systems, however, the question of the optimal choice 
is more complex (see remark in §8). Therefore, we presented here the con- 
struction of the advection scheme with no assumption made about the signs 
of a and [3. 
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5.2. Euler system. The auxiliary variable formulation of the Euler system can 
be written in the following vector form 

s t + U ■ Vs = 0 

(96) pU t + pU - VU + Vp = 0 

Pt + U ■ Vp + pc 2 V • U = 0 


where U = (u, v) 

Consider again triangle T on Fig. 4 and the non-orthogonal coordinate system 
(£,7/) aligning with two of the faces of this triangle (or the unit vectors ei,e* 2 ) Define 
the following quantities 


(97) 
and 

(98) 


U = U-e 2 , V = U-e 1 . 


a = 


U -ni 


P = 


U-n 2 


0 7 ' 0 

We can rewrite the Euler system (96) in this new non-orthogonal coordinates (£, rf) 


(99) 


St + — 0 

plAt + pcdA ^ + p(iU n + p£ = 0 

pV t + paV^ + p(3V, n + Pn = 0 
p t + apt? + (3p v + pc 2 (a^ + fin) = 0 


In order to compute the fluctuation of the system 

(100) 


r = + r v 


on a general triangle T we follow again the two-dimensional linearization procedure 
presented in [16] and assume the quantity varying linearly on the triangle is the pa- 
rameter vector 

(101) m = yfp(l,u,v,H) T 

and its averaged value on the triangle T is given by the following 

~ m 1 +m 2 +m 3 

( 102 ) m = 


Then (similarly to the structured grid case 
(103) 


= —St 


( as£ 

apU ^ + pi 

apV^ 

V apl + c 2 pal j 


where 

(104) 

(105) 

(106) 


Pi 

Pi 


2mi(mi)^ 

7 - 1 


7 


[(m 4 (rai)£ + mi(m 4 )^) + ( m 2 (m 2 )^ + m 3 (m 3 )^)] 


= Pi 


-c pi 
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(107) 


mi(m 3 )t - "^( toi )^ 


P U ( = 


PH 

PH 


r' n can be evaluated in a similar way. Defining the following set of conservative vari- 
ables 


(108) 


u — 


( p \ 

pU 

P v 

V « / 


and denoting R the fluctuations of the Euler system in these variables on triangle T, 
it can be easily verified that the fluctuations in the following relation holds 


(109) 

where 


( 110 ) 


R = C a ■ r 


/ 1 

0 

0 

1/c 2 

\ 

u 

1 

0 

Ufc 2 


V 

0 

1 

V/c 2 


V u 2 / 2 

a 

P 

l/( 7 - 1) + U 2 /( 2c 2 ) 

) 


The fluctuation distribution formulae 


(111) 


s^y 1 

= ul 

+ T 2 C a {rt - rt) 

S 2 u n 2 +1 

= S 2 u n 2 

+ ^C a (r^ + ry 

S 3 u n 3 +1 

= S 3 u n 3 

+ ^C a [(y + r^) + (r v - f 71 )] 


will result in an upwind scheme provided we define the artificial viscosity terms r^,r n . 
However, before doing this we would like to rewrite (111) for the update of the “reg- 
ular” conservative variables u = (p, pu, pv,e) T . Note that 


( 112 ) 


u = 0 • u 


where 


(113) 


/ 1 0 0 0 \ 

0 nf/O n\/6 0 

0 n 2 /0 n^/O 0 

v 0 0 0 1 / 


Introducing the following matrix 


(114) 


/ 1 

0 0 

1/c 2 

\ 

u 

nf/ 6 n\/6 

u/c 2 


V 

n 2 /0 n 2 /0 

v/c 2 


V u 2 / 2 

a (3 

1/ (7 — 1) + U 2 /(2c 2 ) 

' / 


we obtain the fluctuation distribution formulae for the “regular” conservative variables 



Ni< +1 

= Sxul 

+ r_ C o {r i-ry 

(115) 

S 2 u n 2 +1 

= S 2 u% 

+ r_ C o (r v + rV) 


S 3 u n 3 +1 

= S 3 u% 

+ T 2 C e a [(y+ry + (rV-ry] 


16 



The relationship between the artificial viscosity and the fluctuations is given by 

T*^ = M £7*^ 


(116) 


r' n = M r] r' n 


where 

(117) 

and 

(118) 

(119) 


( 120 ) 

and 

( 121 ) 

( 122 ) 


= 


M|“ 6 , if |d| < 7 

, if \a\ > 7, 


Mr, 


/ 

M( ub = 

V 

/ 

= 

V 


I Mf ub , if \(3\ < 7 

: \ Mf up , if |/3| > 7, 

sign(d) 0 0 0 

0001/7 

0 0 sign(d) 0 

0 7 0 0 

sign(/3) 0 0 0 

0 sign(/3) 0 0 

0 0 0 1/7 

0 0 7 0 


\ 

/ 

\ 

/ 


Mp p = sign(d)/, 
M s v up = sign 0)1. 


A linearity preserving scheme can be constructed by introducing the vectors 
r^* ,r v * with their elements defined by 


(123) 


r v* _ r n , i 

i i ' q % ' i 


for i = 1,2, 3, 4, where 
(124) 


ft = 


and substituting them instead of r^,r v respectively in (115), (116). 

Remark 5.1. In order to construct a non- oscillatory linearity preserving scheme 
for the Euler equations, derivatives can be approximated along any two out of three 
faces of the triangle. However, the resolution of such features of the flow as shocks, 
contact discontinuities and shear layers can be affected by the particular choice (see 
§8/. Shocks are being resolved better (i.e. are represented by sharper numerical layers) 
if the numerical derivatives are computed along those two faces of the triangle which 
are the closest to the direction of the shock. The same is true for the shear layers and 
contact discontinuities, although the two faces closest to the discontinuity direction in 
this case will also be either two inflow or outflow faces of the triangle. 
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6. Outline of the extension to three dimensions. A detailed description of 
the three-dimensional schemes both for the scalar advection and Euler system will be 
given elsewhere. In this paper we only outline briefly their construction in the case of 
structured Cartesian grids. Each cube having grid nodes as its vertices can be divided 
into two prisms. The prism can be divided into three tetrahedra each. Any of the 
obtained tetrahedra will have 3 of its edges parallel to x,y and z axes respectively. 
We shall utilize this in this section to make the presentation free of the unnecessary 
details related to the non-orthogonal coordinate systems. The construction of the 
fluctuation-splitting schemes both for the scalar advection and the Euler system will 
be outlined for the case of the tetrahedron T as illustrated on Fig. 5. 

6.1. Advection scheme in three dimensions. Consider a scalar constant co- 
efficient advection equation in three dimensions 

(125) u t + au x + bu y + cu z = 0 

The fluctuation of this equation on the tetrahedron T (see Fig. 5), whose volume is 

V = h 3 / 6. 


(126) 

where 

(127) 


R = R x + R y + R Z 


R x = _/ L _[ a ( M3 _ U4 )] 

Ry = -!^[b(u 2 -u 3 )} 

R z = -irK^i - u 2 )\ 


The fluctuation distribution formulae are given as follows 
(128) 


E < +1 

= Vu^ 

+ %C a (R z + R z ) 


Vu n 2 +1 

= VuV, 

+ %C a [(Rv + R*) + (R Z 

- R z )] 

Vul+' 

= Vu^ 

+ ^C a [(R x + R x ) + (R y 

-R y )\ 

V«J +1 

= Vu'l 

+ T 2 C a (R x - R x ) 



It is easy to see that the following choice of the artificial viscosity 
(129) 


R x = sign (a)R x 
Ry = sign (b)R y 
R z = sign (c)R z 


results in a regular upwind scheme. 

Introducing the following quantities 


(130a) 


(130b) 


(130c) 


Q x = 


Q y = 


Q z = 


mx + m* 

R x + [Ry]t + [R z ]t 

[R x ]~ + [R z ]y 

R y + mt + mt 
m j + m j 

R z + [R x ]f + [Rv]t 
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where 



(131) 

rnt = • 

\ m 
l °. 

and 



(132) 

[VC = I 

0, 

-R\ 

or 



(133) 


m~ K 


if sign (ij >t ) = sign (i? K ), 
otherwise, 


if sign (R L ) = sign (i? K ), 
otherwise. 


= mi - r 1 


where i and k stand for one of (though different) x,y or z each. 

Then we define the following quantities 

R x * = R x + V(Q y )R y + 'S(Q Z )R Z 
(134) R y * = R y + ^>{Q X )R X + V{Q Z )R Z 

r z * = r z + ^(g^)^ + 

where T is a non-compressive limiter. 

Substituting R x * ,R y * ,R Z * into (128), (129) instead of R x ,R y ,R z respectively we 
obtain a linearity preserving (second order accurate in this case) positive advection 
scheme. 

To demonstrate this we can assume without loss of generality that 


(135) 

Q y -Q z > 0 
Q x (Q y + Q z ) < 0 

Then 


(136) 

Qy = -R x /(Ry + R z ) 
Q x = l !Q y 
Q z = Q y 

In this case 


(137) 

^(Q y )R y + y(Q z )R z = 
^(Q X )R X = -^(Q y )(R y + R z ) 

and therefore 


(138) 

R x * = (1 -V(±))R X 
R y * = (1 - ^{Q y ))R y 
R z * = (1 - ^{Q Z ))R Z 


It can be concluded that if T is a non-compressive limiter, then the constructed scheme 
is positive. 

If the fluctuation R vanishes then 


(139) R x = R y + R z 
and, therefore, 

(140) Q x = Q y = Q z = 1. 

Then it is easy to see that 

(141) R x * = R y * = R z * = 0, 

provided T(l) = 1, i.e. the constructed scheme is linearity preserving. 
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6.2. Euler system. Euler equations in three dimensions 
(142) u t + F(u) x + G(u) y + H(u) z = o 

where 


u = 


( p \ 

pu 

pv ; 
pw 

\ e 


( P u \ 


( P v \ 


( p w \ 

pu + p 


puv 


pwu 

puv 

; G(u) = 

pv 2 + p 

; H(u) = 

pwv 

puw 


pwv 


pw 2 + p 

\ u(e + p) y 


\ v(e + p) y 


\ w(e + p) j 


Introducing the auxiliary variables (s, u, v, w,p ) we can rewrite these equations in the 
following non- conservative form 


S t + us x + VSy + WS Z = 0 
put + puu x + pvUy + pwu z + Px = 0 
(143) pv t + puv x + pvVy + pwv z + Py = 0 

pw t + puw x + pvWy + pww z +p z = o 

Pt + UPx + VPy + wp z + pc 2 (u x + Vy + U) Z ) = 0 


Fluctuation on tetrahedron T (see Fig. 5) is 

R = J J J u t = - J J J {F x + G y + H z )dx dy dz 

(144) 


= -V T [F X + G y + H z 


Fluctuation in auxiliary variables on each tetrahedron can be represented as a sum of 
three parts 


(145) 


r = r x + t v + t z 


The matrix C a relating the fluctuations in conservative variables R to those in the 
auxiliary variables 


(146) 


C a 


/ 1 0 0 0 

u 10 0 

v 0 10 

w 0 0 1 

y tj 2 / 2 u v w 1/(7 


1 /d 2 \ 

u/c 2 
v/c 2 
w/c 2 

1 ) + U 2 /(2~c) 


It was pointed out in [16] that the two-dimensional conservative linearization presented 
there extends to three dimensions in a straightforward way. Again, we assume that 
the quantity that varies linearly on the tetrahedron T is the “parameter- vector” and 
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the averaged quantities can be defined in a manner similar to the two-dimensional 
case. The fluctuation in the auxiliary variables on T can be written as follows 


(147) r = -S T ■ A ■ r x = -S T A ■ 

The fluctuation distribution formulae in this case are 


(148) 


Similarly to the two-dimensional case, the relationship between the the artificial vis- 
cosity needed in ^-direction to obtain the upwind scheme and r x is the following 


Vu" +1 

= Vu™ 

r 2 Ca{r z + r z ) 



Vu n 2 +1 

= Vu n 2 

+ r'a[(r y +r y ) 

+ 

(r z — f z )\ 

Vu n 3 +1 

= Vu% 

+ ^C a [{T x + f x ) 

1 + 

( r y — f y )\ 

Vu n 4 +1 

= Vu'l 

+ T 2 Ca(r x ~r x ) 




(149) 


r x = sign(A)r :r 


r y and r z can be evaluated in a similar way. 
Denoting 


M x = sign(i) 


it can be verified that 

(150) 

where 


M x = 


M s x u \ 

Mr, 


if \u\ < c 
if \u\ > c, 


(151) 

and 

(152) 


M sub = 


( sign(fi) 0 

0 0 
0 
0 

\ 0 c 


0 
0 

0 sign(fi) 

0 0 
0 


0 

0 

0 

sign(fi) 

0 


0 \ 

1/7 

0 
0 

0 / 


M s r = sign (u)I 


Following the procedure of constructing a high resolution scalar advection scheme in 
three dimensions, presented briefly in §6.1 we define the following quantities 


(153) 


% = 


£ 

rf + [r y t ]t + [r z ]t 


(154) 


y 

9i = 


[rf]y + [r:. 


i \y 


r y + [rf]t + [rf]t 


(155) 


Qi = 


[rf] z + [r. 


y }~ 
J Z 


rf + [rf]+ + [r y ]f 
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where [. . .]±, [. . ,]±, [. . ,]f are defined by (131), (132). Now we can introduce 

rf = r f + '&(<li) r i + ^(^i) r i 

(156) r\ = r\ + ^(qf)rf + T(gf)rf 

r z* = r z + q( q f) r x + T(gf)rf 

for i = 1, 5. 

Substituting r x * ,r y * ,r z * instead of r x ,r y ,r z into the fluctuation distribution for- 
mulae (148) and (149) we obtain a linearity preserving (second order accurate for the 
case of structured meshes) genuinely three-dimensional scheme for the Euler equations. 

Note that 

(157) r = r x + r y +r z = r x * + r y * + E* . 

Therefore, the constructed scheme is conservative. It is also easy to see that it is 
linearity preserving. 

7. Multigrid solver. One of the most attractive properties of the constructed 
multidimensional scheme for the Euler equations is that the Collective Gauss- Seidel 
relaxation is stable when applied directly to the high-resolution discrete equations. 
This can be utilized to construct a very simple and efficient steady-state multigrid 
solver. In this work we would like just to illustrate the basic advantages of the multigrid 
solver that uses the constructed multidimensional scheme. Therefore, we do not discuss 
here any issues concerning the multigrid algorithm for general unstructured meshes, 
but address the interested reader to the literature (see, for instance, [11]). 

We shall present briefly in this section the simple multigrid algorithm used in the 
numerical experiments presented in this work. 

Description of the algorithm. The multigrid cycle used in this work is IT-cycle 
and the entire algorithm can be implemented either as Cycling or as Full Multigrid 
(. FMG ). 

Relaxation. The algorithm employs the Collective Gauss- Seidel relaxation as a 
smoother. In order to update the solution at each node a system of four nonlinear 
equations has to be solved. This can be done using the Newton method. One iteration 
at each node is sufficient for this purpose (see [18] for the discussion on this matter). 
The ordering of the relaxation can be Red-Black, lexicographic etc. 

Restriction and prolongation. Assume that we have a hierarchical sequence of tri- 
angular grids, i.e. that any triangle on the coarser grid is a union of four triangles 
on the finer grid. The natural choice for prolongation in this case would be a linear 
interpolation along each face of the triangle belonging to the coarser grid. The re- 
striction in this case can be just combining the fluctuations from the four finer grid 
triangles forming the coarse grid one. However, we consider only Cartesian grids in 
the numerical experiments presented in this work (see Figs. 2, 3). In this case we can 
use more conventional prolongation and restriction operators: bilinear correction in- 
terpolation and Full- Weighting of the residuals (see [1]). Term “residual” is used here 
in its usual sense: residual of the discrete equation at each gridpoint (constructed of 
the contributions from the triangles having this gridpoint as a common vertex). 
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Previous relevant results. Application of multigrid methods for the advection equa- 
tion was studied in details in [17], [18] in conjunction with the genuinely two-dimensional 
control- volume type scheme developed there. It was demonstrated there that applying 
2FMG — W(2, 1) algorithm is sufficient to obtain high quality solutions (i.e. solutions 
with sharply resolved discontinuities and second order accuracy both in smooth regions 
and in discontinuity location) for the advection equation. Many of the conclusions of 
[17], [18] apply directly to the Euler system case considered in this work. 

Efficiency of the multigrid solver for advection-dominated problems. It was shown 
in [1] that the residual reduction per multigrid cycle for an advection (dominated) 
problem cannot in general be better than .75 for the second order accurate discretiza- 
tion (.5 for the first order). This is because the coarse grid provides only a fraction (.75 
or .5 respectively) of the needed correction for some components. The ability of the 
multigrid algorithm to achieve rate of convergence close to .75 (for the second order 
accurate approximation) means that the obstacle towards achieving a better efficiency 
are not the smoothing properties of the relaxation but rather an unsufficient coarse 
grid correction for certain components, (see §9.2.2 for further discussion on this issue). 

8. Numerical experiments. The purpose of the numerical experiments re- 
ported in this section is to verify the robustness of the constructed scheme and the 
quality of the numerical solutions obtained by its means. Some preliminary experi- 
ments illustrating the performance of the multigrid algorithm that utilizes this scheme 
are presented as well. 

8.1. Supersonic flow in a channel with a bump. The test case considered 
here is a supersonic (Mach=2.9) flow in channel with a circular bump. The bump is 
located at the lower wall of the channel at 1 < x < 2. and its surface is a circular arch 
of 7t/3 and radius 1. Note that the actual shape of the domain is a rectangle. The 
influence of the bump on the flow is imposed through the boundary conditions: the 
velocity component normal to the surface of the bump at a certain location is being 
reflected. 

The first experiment uses a grid of the size 200 x 40 points. The density contour 
plots of the steady-state solution are presented on Fig. 6(a). The scheme used is the 
one given by (71), (72), (73) in §4 with the minmod limiter. 

The second experiment presented on Fig. 6(b) corresponds to the same settings, 
except that the grid is twice finer (400 x 80 points). As it is expected, the grid 
refinement results in a better resolution of the flow features. 

The third experiment (Fig. 6(c)) is performed on the grid of the same size as 
the first one. However, the triangulation is as illustrated on Fig. 3 and the scheme 
employed is constructed according to (115), (116), (123) in §5 and the derivatives are 
approximated using the diagonal and parallel to the x direction faces of the triangles. 
Even though the grid is structured, the most essential feature of the unstructured grid 
formulation of the scheme - use of the non-orthogonal coordinate system is present. 

The purpose of this experiment is to test the performance of such a scheme and to 
illustrate the effect of the scheme and triangulation on the resolution of discontinuities. 
It can be seen on Fig. 6(c) that the shock reflected from the upper wall is resolved better 
than the stronger one that is incident to the upper wall. This is because for this 
scheme/grid combination discontinuities whose direction is close to the grid diagonal 
(upper-left to lower-right) will be resolved very well (in addition to those which are 
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close to the horizontal or vertical directions). 

8.2. Transonic flow over a circular bump. The testcase considered here is 
a transonic flow (free-stream Mach= .9) over a flat wall with a bump (Fig. 7). The 
surface of the bump is a circular arch of 7t/3 and radius 1 and its location is between 
3.5 < x < 4.5. Again, in order to keep the experiments simple at this stage of work, 
the bump is treated the same way as in the previous experiments. The grid is 200 x 200 
points and the scheme used is presented in §4. The shock of the “fish-tail” shape can 
be clearly observed on Fig. 7. 

8.3. Low Mach number flow over a circular bump. Here we present a 
numerical experiment concerning a low Mach number ( = .l) flow over a flat wall with 
a circular (arch of 7t/3 and radius 2) bump. Here as well as in the previous case 
the presence of the bump is imitated through the appropriate boundary conditions. 
The grid is 200 x 200 points. The density contours of the steady-state solution are 
presented on Fig. 8. 

8.4. Multigrid algorithm. To illustrate the performance of the multigrid algo- 
rithm we consider here the well known testcase of a shock reflecting from a flat wall. 
The multigrid algorithm involves five grids (levels): the finest consists of 129 x 33 
points, the coarsest -9x3 points. 

The multigrid algorithm is based on the scheme presented in §4 used with the 
lexicographic Gauss-Seidel relaxation. The restriction and prolongation procedures 
are the standard Full Weighting of the residuals and bilinear correction interpolation. 
The first experiment of the first concerns performing one 1T(2,1) cycle. Fig. 9(a) 
presents the initial guess and Fig. 9(b) - the numerical solution obtained by one W(2, 1) 
cycle. It is remarkable that such a little computational work is sufficient to make the 
reflected shock is clearly visible. The numerical solution to this problem obtained by 
2 FMG — W(2, 1) algorithm is presented on Fig. 9(c). 

Note, that in this case the flow is aligned with the ^-direction in the significant 
part of the domain. In this case the artificial viscosity in the cross-stream direction in 
the entropy and n-momentum equations (see §4) vanishes. Therefore, no smoothing 
can be obtained in y-direction in some components. A multigrid algorithm utilizing 
the time-stepping type relaxation can deal with such a situation only using the semi- 
coarsening technique. Our algorithm employs the Gauss-Seidel relaxation. Therefore, 
it offers a much simpler and more efficient treatment of this problem: relaxation with 
lexicographic ordering in the stream direction. 

The rate of convergence observed in this testcase as well as in other simple ex- 
periments concerning variety of flow regimes is very close to .75 (see §§7,9. 2. 2 for a 
discussion). 

9. Discussion and conclusions. 

9.1. Summary of the current work. A new two-dimensional high-resolution 
(at the steady-state) scheme for the compressible Euler equations was presented. It 
is triangle-based and can be formulated with the same degree of simplicity both on 
structured and unstructured grids. The main advantage of this scheme is that Gauss- 
Seidel relaxation can be applied directly to the resulting discrete equations. This 
allows to construct a simple and efficient multigrid steady-state solver. 

A remarkable property of the constructed scheme is also its very compact stencil: 
it involves only the immediate neighbors of the point of interest. 
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A variety of the flow regimes (supersonic, transonic and low Mach number flow) 
were considered in the numerical experiments to verify the quality of the solutions 
obtained by means of the constructed scheme and to demonstrate the efficiency of the 
multigrid algorithm. 

Generalization of this scheme to three dimensional tetrahedral meshes is presented 
briefly as well. 

9.2. Future work. 

9.2.1. Compressible Navier-Stokes equations. The extension of the con- 
structed scheme to the case of compressible Navier-Stokes equations is straightforward. 
It is interesting to note that the artificial viscosity present in the momentum equations 
in the subsonic case contains already terms of the type (u x + v y ) x and (u x + v y ) y that 
appear in the physical viscosity of the compressible Navier-Stokes equations. 

9.2.2. Further improvement of the multigrid efficiency. As it was men- 
tioned in §7 the main obstacle towards the further improvement of the multigrid 
efficiency is the following fact: for the hyperbolic problems the coarse grid correction 
is not sufficient for certain error components. 

This difficulty was already addressed in the literature and some techniques to 
improve the multigrid efficiency were developed in [2]. Therefore, one possibility is to 
adapt these techniques for our case - compressible Euler equations. 

We shall mention here another way to deal with this difficulty. The steady-state 
compressible Euler equations in subsonic case are very similar from the mathematical 
point point of view to the incompressible Euler equations - both problems are of the 
mixed elliptic-hyperbolic type. 

There exist very efficient steady-state multigrid solvers for the incompressible flow 
equations, which rely on the primitive variables formulation (velocities and pressure). 
Most of them employ staggered grid discretization of the equations. A typical efficiency 
of such algorithms for a low Reynolds number case is comparable to that for the Poisson 
equation. Although the converegence rate deteriorates for the advection dominated 
(high Reynolds number) case, it may still be possible to recover the optimal efficiency 
by treating the advection factor in some special way (like relaxing in the flow direction). 
This is achieved due to a certain “separation” of the elliptic and hyperbolic factors of 
the equations in the treatment (see [1]). 

However, extending this approach to the compressible Euler equations appeared 
to be a very difficult problem despite the similarity of the equations. 

A canonical variable formulation of the Euler equations for both incompressible 
and compressible cases was suggested recently by Ta’asan [23]. This formulation 
facilitated the construction of the discrete scheme (using staggered grids) and of the 
relaxation procedure that separate a treatment of the advection and full-potential 
factors (see [24]). The resulting multigrid algorithm achieves a very good efficiency in 
subsonic case if the relaxation is performed in the flow direction. 

One of the future directions is to modify the scheme presented in this work so that 
it will allow to achieve a similar efficiency for the case of triangular (non-staggered 
and, possibly, unstructured) grids. This may be possible because some very impor- 
tant similarities between the present scheme and the existing discretizations for the 
incompressible flow equations used in very efficient multigrid solvers (see, for instance, 
[1]) can be observed (see §9.3). The resulting discretization will rely on the primitive 
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rather then the canonical variables. The latter is crucial for the further generalization 
of the solver to the compressible Navier- Stokes equations. 

9.2.3. Time-dependent problems. The constructed scheme is capable of pro- 
ducing high-resolution (second order accurate on the structured grids) steady-state 
solutions for the compressible Euler equations. A very natural way to extend this 
capability to the transient problems is to use a Lax-Wendroff type modification of the 
constructed scheme, i.e. to scale the artificial viscosity in such a way that it will cancel 
the time error of the forward Euler time-stepping due to the equation (or system of 
equations) being solved. 

To illustrate this on the simple case of the two-dimensional advection consider the 
scheme for the structured grids presented in §2.1.2. The use of the following artificial 
viscosity 

R x = jaR x 

Ry = £ bRy . 

instead of one given by (5) results in the high-resolution scheme (constructed by sub- 
stituting R x * ,R y * defined by (6) instead of R x ,R y ) that is second order accurate in 
time and space. However, it is easy to see that in in general this scheme is not of the 
positive type anymore (see also [6]). 

For the case of the Euler system modifying the construction of the high-resolution 
scheme presented in §4 by substituting the artificial viscosity defined as follows 

r x = lAr x 
r y = IBry 

instead of (72), results in a Lax-Wendroff type scheme for the Euler system. However, 
the applicability of this scheme is restricted to the problems with with no strong 
shocks. 

It is interesting to note that the two-dimensional schemes by Colella [3], [4], LeV- 
eque [10], Radvogin [12] can be interpreted as the schemes of the Lax-Wendroff type. 
However, the nonlinear mechanism that enables them to deal with (strong) shocks 
relies on one- dimensional limiters and characteristic variables in x and y direction, 
which introduces an essential flavor of dimensional splitting. 

The first step towards the construction of a genuinely multidimensional time- 
space accurate and robust Euler scheme should be to construct a positive time-space 
accurate scalar advection scheme of the type Lax-Wendroff type. This is one of the 
direction being currently investigated. 

9.3. What is a truly multidimensional scheme ?. Perhaps, it will be easier 
to answer this question after the nonoscillatory time-space accurate multidimensional 
scheme for the compressible Euler equations is constructed. However, some attributes 
of the truly multidimensional schemes can be seen already on the example of the 
scheme constructed here. 

One of the main characteristics of a truly multidimensional scheme for fluid dy- 
namics is that the divergence operator (in the pressure equation, the auxiliary variables 
formulation of the Euler system (56) should be represented in a difference scheme as 
whole. This is not the case for the dimensionally split schemes. It can be clearly seen 
on the example of the dimensional upwinding scheme (71), (72) that terms u x and 


(159) 
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Vy are separated (when contributing to the artificial viscosity of the momentum and 
pressure equations). 

Consider the case of subsonic flow. In this case the fluctuation of the pressure 
equation contributes to the artificial diffusion of the momentum equations. Let us 
modify slightly the scheme constructed in §4. and compute the fluctuation of the pres- 
sure equation contributing to the artificial viscosity terms of the momentum equation 
according to the following 

(160) rf = r\ = r 4 + rf 

The only difference between (160) and (73) is that the limiter function is set to unity. 
Although 


rf + rf = 2 r 4 ^ r 4 , 

such a scheme is conservative, since (160) is used only in the artificial viscosity. This 
is because the role of the artificial viscosity can be interpreted as “redistribution” of 
the fluctuations split by the central scheme - to subtract a certain amount from one 
vertex of the triangle and to add the same amount to another vertex. This does not 
violate conservation, regardless of what amount has been “redistributed”. 

Recall, that for some very popular schemes for incompressible flow computations 
(see, for instance, [1]) the residual of the continuity equation contributes to the veloc- 
ities update in such a way that the discrete vorticity remains unchanged. 

The “physical” meaning of this is that the change in (evolution of) the vorticity 
is only due to the residuals of the momentum equations, as one would expect. 

“Mathematical” meaning of this fact is that the elliptic factor of the incompressible 
Euler equations is being discretized and relaxed as such. This allowed to construct a 
highly efficient multigrid solver (see [1]). 

What can be observed clearly in the subsonic case is that the fluctuation of the 
pressure equation (being a part of artificial diffusion of the discretized momentum 
equations) contributes to the change of the velocities (divergence update) according 
to exactly the same pattern as the residual of the continuity equation in the schemes 
for discretizing the equations of incompressible flow. 

Thus, the scheme presented in this work establishes a link between some standard 
techniques used for solving equations of the incompressible flow and the class of upwind 
schemes commonly used for compressible flow computations. 
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Fig. 6. Supersonic flow in a channel over a circular bump: a) grid 200 x 40 pts., triangulation I, 
scheme from §4; b) the same, except the grid 400 x 80 pts.; c) grid 200 x 40 pts., triangulation II, 
scheme according to §6, which relies on the approximate derivatives along the nonorthogonal faces oj 
the triangles. 





Fig. 8. Low speed flow (Mach= .1) over a wall with a circular bump. 
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